{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000000/1000000 [00:02<00:00, 365191.56it/s]\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from tqdm import tqdm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "NUMBER_OF_STEPS = 1000000\n",
    "SAMPLE_INTERVAL = 1\n",
    "\n",
    "length = int(NUMBER_OF_STEPS/SAMPLE_INTERVAL)\n",
    "x_list = np.zeros(length, dtype=float)\n",
    "v_list = np.zeros(length, dtype=float)\n",
    "inv_list = np.zeros(length, dtype=float)\n",
    "vel_eta_list = np.zeros(length, dtype=float)\n",
    "\n",
    "def nh(pos_eta, vel_eta, mas_eta, Ek2, kT, dN, dt2):\n",
    "    dt4 = dt2 * 0.5\n",
    "\n",
    "    G = Ek2 - dN * kT\n",
    "    vel_eta[0] = vel_eta[0] + dt4 * G\n",
    "\n",
    "    pos_eta[0] += dt2 * vel_eta[0] / mas_eta\n",
    "\n",
    "    # Compute the scale factor\n",
    "    factor = math.exp(-dt2 * vel_eta[0] / mas_eta)\n",
    "\n",
    "    # Update thermostat velocities from 0 to M - 2:\n",
    "    G = Ek2 * factor * factor - dN * kT\n",
    "    vel_eta[0] = vel_eta[0] + dt4 * G\n",
    "\n",
    "    return factor\n",
    "\n",
    "\n",
    "# Main simulation loop\n",
    "dN = 1.0\n",
    "dt = 0.01\n",
    "dt2 = dt * 0.5\n",
    "kT = 1\n",
    "x = 0\n",
    "v = 1\n",
    "mass = 1.0\n",
    "k_spring = 1\n",
    "f = -k_spring * x\n",
    "pos_eta = [0.0]\n",
    "vel_eta = [1]\n",
    "mas_eta = 1\n",
    "\n",
    "Ek2 = 0.0\n",
    "factor = 0.0\n",
    "\n",
    "for step in tqdm(range(NUMBER_OF_STEPS)):\n",
    "    inv = mass * v * v * 0.5 + k_spring * x * x * 0.5\n",
    "    inv += kT * dN * pos_eta[0]\n",
    "    Ek2 = v * v * mass\n",
    "    factor = nh(pos_eta, vel_eta, mas_eta, Ek2, kT, dN, dt2)\n",
    "    v *= factor\n",
    "\n",
    "    v += dt2 * (f / mass)\n",
    "    x += dt * v\n",
    "    f = -k_spring * x\n",
    "    v += dt2 * (f / mass)\n",
    "\n",
    "    if step % SAMPLE_INTERVAL == 0:\n",
    "        nn = int(step/SAMPLE_INTERVAL)\n",
    "        x_list[nn] = x\n",
    "        v_list[nn] = v\n",
    "        inv_list[nn] = inv\n",
    "        vel_eta_list[nn] = vel_eta[0]\n",
    "\n",
    "    Ek2 = v * v * mass\n",
    "    factor = nh(pos_eta, vel_eta, mas_eta, Ek2, kT, dN, dt2)\n",
    "    v *= factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAFpCAYAAABu98hvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAngUlEQVR4nO3dfXRU1b038O+PQHwBbUFCeQqUsKKXFN8oRvClpbTSXl685tKH9iG3xVJLs8CX2mpXDbW9tdVbcbVVH1qES1Gp1SuPbeGKEktNW8TeKiXQoCKJiyjIi5ZRLJZgCSG/549MhpnJJDkz2Wf2Pnu+n7VYK+fk5Jxfwi/fnNlzzj6iqiAiomjrZ7sAIiLqO4Y5EZEHGOZERB5gmBMReYBhTkTkAYY5EZEHjIa5iBSJyF9E5EmT+yWyjb1NrjN9Zn4jgJ2G90nkAvY2Oc1YmIvISAAzAaw0tU8iF7C3KQpMnpnfC+CbANoN7pPIBfeCvU2O629iJyJyJYCDqrpVRKb0sF01gGoAGDhw4EXl5eUmDk+ErVu3vqWqJab3y94m24L2tpiYm0VE7gQwF0AbgFMBnAlgjap+obuvqaio0Pr6+j4fmwgARGSrqlaEsF/2NlkVtLeNDLOo6iJVHamqpQDmAPh9T81OFBXsbYoKXmdOROQBI2PmyVR1I4CNpvdLZBt7m1zGM3MiIg8wzImIPMAwJyLyAMOciMgDDHMiIg8wzImIPMAwJyLyAMOciMgDDHMiIg8wzImIPMAwJyLyAMOciMgDDHMiIg8wzImIPMAwJyLyAMOciMgDDHMiIg8wzImIPMAwJyLygPFngFLuSmvWZ7X97sUzQ6qEyKyNjQdx/aPbcPxEO461KQYUAf37CYr6CVqOtWPI6QPwt38cxx2V56Fq0mjb5UYSw9yysYvW45jm9rWd4T9wALDjdgY7uWdJXRPurtvVZf3xE8DxEwqgo/nfPnocALBo7Uv41tqX8ImxJfjR58ZjyMDifJYbaQxzS7I9C+9Jy/GO/fUHsItn6+SAbXvewYJf1OPgkdaU9UX9gBPtyHhm3hnoCuD3TTFc+oM6/OfcCkwpH2bhO4gehnme3VX7MpZtei2UfbehI9TvnMWXqmRPc+wIPr/yebx3vD2xrrgf8K8TRqJm+oe7Pdtujh1B9UP1aI61AACOnVDc9MsGbPvOp/NSd9SJao6v8fuooqJC6+vrrRzbFpNn40EU0pi6iGxV1QrbdQCF2dtAx9n4/J9vwaH4GfaAfkC/fv3wo9kX4F/Gjwi8n42NB3Hdf21FPxH8YNb5WP/iG6jb+VfcXqDj6UF7m2fmeZLvIO88ZiEFOtlV/dDJIAeA33z94ygrGZT1fqaUD8OO708HAPznM834zY6/AugYT2851ob5k8vMFOwZXpqYBzaC3IVjU2HY2HgQ5333N3ir5WSQf3tGeU5Bnu6zFaNw+oCixPIdtY14dPOePu/XRwzzEFUuecaJMC2tWY+Vm5ptl0EeeqJhP+at2oIjx04k1t056zxjZ89DBhbj4fmTMGTggMS6RWtfYqBnwDAPyeylz2L7gSO2y0i4o7YRTzTst10GeWTbnndww+qGlHWr5l1sfFx7wujB2PadT+P0AZJYt2jtS9jYeNDocaLOWJiLyCgR+YOI7BSRHSJyo6l9R838Bzejfu+7tsvo4obVDfwFyBL7OrNDLa2Ye//zKet+Mmd8qJcR3vf51PcA563agm173gnteFFj8sy8DcDNqvphAJcAuE5ExhncfyQ80bAfdU1v2S6jW/NWbbFdQtSwrzP4xmMNaGk9eenhqnkXZ3XFSi6mlA/DmoWXpayb98CfQz1mlBgLc1V9Q1W3xT/+O4CdAML933VQ+stOF7kwjh8V7OuultQ14fdNscRy2GfkySaMHoybpp6dWH73WBtqftWQl2O7LpQxcxEpBfARAJvD2L+rohSSUarVFYXa18ke3bwn5fb8T44tCf2MPN1Xp47FvEs/lFheXb8fS+qa8lqDi4yHuYgMAvBrAF9T1XfTPlctIvUiUh+LxTLvIKKiGI5RrNmWnvo6/nlve7vToZZWLFr7UmK5SIAffW68lVpuqzwfHzjj5J2kd9ftKvj3g4yGuYgMQEfDP6Kqa9I/r6orVLVCVStKSkpMHtqqKIfix+6ss12C83rra8Df3k72pQdTX5Dc/8WLrU6E9V/Vl6Ysf+Whwn4/yOTVLALgfgA7VfVuU/t1XdRf3u09fMx2CU4r1L5Od1fty9i+7+QLkoWTx1ifAKusZBDunHVeYvl4Owp6/NzkmfnlAOYC+KSINMT/zTC4fydlmt4zaqL8yiIPCrKvk21sPJgyOdyg4iLcMsONC3qqJo3GhR88I7G8un5/wQ63mLya5Y+qKqp6gaqOj/+rNbV/F/kUgj59LyYVYl+nm//z1OGLh748yVIlmT345UtSluet2oLmmDs37OUL7wDN0UXf+43tEohC9/VHt6ItaWLVm6aejQmjB9srKIMhA4vx7RnlKes+t+x/LFVjD8M8R2+/d6L3jSKGZ+eUbEldE9ZufzOxPOvC4fjq1LEWK+re/MllmHXh8MTy20fbIv9+VrYY5jnwOfR8/t4ouObYkZT3g/r3A+6pushiRb1Lr+/uul041NLazdb+YZhn6a7al22XQBS6zyx9NmV55dUXW6okO+nDLXOWF85wC8M8S2E98s0lPDsvbF9/dCsO/+PkvCuzLhxu/TLEoNKHW16JHcXXH91qsaL8YZhngSFHvtvYeDBlnBxwf3glXXq9a7e/WRCXKzLMKSP+4SpM16TNqvmTOePtFNJHyTcTAcD8Arg7lGEeUCGG2/wHC3Y+qYJ09c+eQ3vS8rxLP5T3SbRMqZo0GnMqTtbe1g7vh1sY5tQtl+dlJ7Nue/xFbGo+lFg+89Qi3FZ5vsWK+m7x7PEpj5tbu/1Nry9gYJgHUIhn5Z04EZf/mmNHsOq511PWrb3uo5aqMeuXC1IfZuHzBQwMc+oRJ+Ly32fuS71876apZ6OsZJClasxKn4wL6Hg+r48Y5r0o5LPyTrc9/qLtEigkS+qacPi9tsTy1LFDnb3LM1dVk0aj9KzTEsv1e9/Fo5v3WKwoHAxz6lX6S3DyR/qsnyu/5NYkWqasuTZ12Cj5IRu+YJj3gGfl5LP0ub+jehliEEMGFmPh5DEp667+2XOWqgkHw5wC4R82v2xsPIjV9fsTy+M+MDCylyEGdcuMcRhULInlTc2HvJoql2HejaoCmtOBCk/6TTQPV1/WzZZ+efyGySnL6XPQRBnDvBvP7f6b7RKc8088O/fCkromtCXdHTSnYoTVZ3nmU1nJIMw49+Q8M4f/0e7N2Xl/2wVQdBTOZKJ+S3/Tc/Hs8Xk5bpChut2LZ4Zex31zL06pZfZ9f8Rfvjst9OOGjWGeAceHyVfpE06lTxkbhmx+nzq3DTvUZ5w7DLU7On4W73jyoBkOs1BW+Icu2r6SNlY+f3JZaMeaec/GnPultGY9rv1FeJNj3Tc3dX52H+Zt4Zm54zrPUBiiZMLxpLHy5Hm/TfvnH/8BTbGjGT+Xfta9sfEg5q3qGty1Ow7i2l9s6RK8plzwwTPwwoG/A+iYt+WeqlAOkzcMc0elNzxDnfoq/a7HsOYpP9TSmjHIf3fzxzNOEzClfBh2L56J5tgRXPHjZ1I+1zkUEoZVX74EE25/OrT95xuHWdK4EJY9jRfm4w0i8tO38nTXY6aA3L14Zq/zvZSVDMrY32H9TqZfwRP1oRaGuWOChLXtQHfhDx5lT5M+Tp98ypRMT/TJtl/z2d/Jc56nP2EpahjmDsmmiW0HOkVb1aTRoew3fez7dzd/PKf9pP+xCesEIl+XZeYDw9wRDGcKk61ZAnOdSjesPzY+Y5hHGP8AUFA2Zgnsa3/amPgrffKxKGGYJ7H1SKm+NP25HxhosBIqBGGNl5seCkmf+KtyyTPdbNk3pWednvg4efKxqDEW5iIyTUSaRGSXiNSY2m8+RfGRUuu/PsV2Cd7zobeTRXUIY/uBcOZQuX9eONex55uRMBeRIgBLAUwHMA5AlYiMM7Fv33GoxG3s7dzcNPVsI/sZfGr4gwe+PCLP1E9qIoBdqvqqqrYCWA2g0tC+qRf8gxAq9nYOTD167i+3TTeyn0JgKsxHANibtLwvvo4o6tjbFAmmwlwyrNMuG4lUi0i9iNTHYjFDhyYKFXubIsFUmO8DMCppeSSAA+kbqeoKVa1Q1YqSkhJDhyYKFXubIsFUmG8BcI6IjBGRYgBzAKwztG8im9jbFi2pa7JdQmQYCXNVbQNwPYANAHYCeExVd5jYN/WOc6WEh72dm8v+w8xshOlPRaLuGbvuR1VrVfWfVLVMVf/D1H59xyB2n2+9fagl/AcAHvi7+WOEdZGirakOTOMdoEmieIkf/xhQtuav+nMo+x31vlNC2W+nV0P6/Uye6iDTu91RwTB3AAOZwnbGqUWJj7ftPRzKMZ5dNDVlua+Pffvwrfn/vfhBSFMd5APD3BG5BDr/CFBQP//SpLwfs69PCbLxnOWoTnUAMMwji0FO2ZgwenBejlOUtnzx9zfktJ/0/g5rQrlte94JZb82MMwdEjSgbQd5FN9boFTzH9wcyn6b03ojdrQt68sLP3ZnXZd1YU0oN3fl86Hs1waGuWN6C2rbQU7RNevC4YmP65reCu04xWnLd9ftwtU/ey7Q1178/Q3Ye/hYyrqKUWcaqqyrluPtiY/nXfqh0I6TDwzzNC6cdZbWrM8Y2gxy6ot7qi5KWQ5riOGVDL9Dm5oPBTpRiR1t67L+V9d9zFhtyVZuak5Zvq3y/FCOky/9bRdA3WN4U5jmPfBnvPC9fw5l37sXz+zxhEQAvNbNNun7CcsdtY2h7dsGnplTVlx45UK5+/aM8sTH7x5rQ3MsnAc+AD33iqLnkxXp5ev7amNj6pU2YT19KZ8Y5hkwsMhX8yeXpSx/dvmfQj3e7sUzMXXs0Ky+5tszyvFayL+D81elXgMf5UsSO3GYhajAzLpwONZufxMAcKjlOA61tGLIwPS3Lc1ZmXSN+22Pv4hVz73eZZtvzyjv8ocmLEvqmpA8Mj+nwo/p6RnmFBhfsfjhnqqLsHb7ySGO6p9vwa+uvTwvx76t8nzrbzSmT961ePZ4O4UYxmGWbjC4yGfJlynWv/63Lld2+Oqu2pdTlk09q9QFDHMKJMoTEFFX6Zcp+nZlR3eWbXotZdnUs0pdwDDvwdiS022X4Iyw35Ci/Es/K73t8RctVZIfs5c+m7LswxUsyRjmPdhw8ydsl0AUmq9OHYvTB5xcXvXc694+2eeu2pdRv/fdxPLksiFeXMGSjGFOveL7B/56eP5lKcs+Ptln2553ugyvPPSVSy1VEx6GeS8YZOSzCaMHY+HkMSnrgs6jEhWfT/t+fBte6cQwpx7xj5n/bpkxDmcljbdsaj7kzdUtNb9qwHttmlj2cXilE8M8AAYa+e6xhanDLXfUNublWaFheqJhP1bX709Z5+PwSieGOXWLf8QKR1nJoJR5WwCgavn/WKrGjBtWN6Qs+zq80olhHhCDjXw3f3IZSgefmlhuih2N7NUt6c8fnVMxwtvhlU4M8ywk3zXnO/7xKkxrrk+dO/zuul1dZhh03aOb96Q8f/Ss0/t7c8t+TxjmWUi/a85Xq+ZdbLsEsmTIwOIuwxHzf76lm63dc6ilFYvWvpSy7rGF+Zl3xjaGeZYK4Yx1Svkw2yWQRVWTRqfMJNimXec0cdXclamXIS6cPAZlJYMsVZNfDHNKUQh/rKh3i2ePR/+kCXmWbXrN+eGWlZuaseONkw/bmHHuMNwyY5zFivKLYZ4DXwPvjPCmtKYIWvnF1OG26ofcHW5pjh1JmSysfz/gvrmFNVzIMM/Rtu98ynYJxr34fT//SFFuppQPSxluaW13dzKuzy77Y8ryyqsLK8gBQ2EuIj8UkUYReUFE1orI+03s12VhPpnFBl9fbfRFIfZ1usWzx2P4mackllc99zq27XnHYkVd1fyqAYeOnkgsL5w8piDf9zF1Zv40gPNU9QIArwBYZGi/TvMlAC8tfb/tElxVkH2d7pGvXJKy/IWfPW+pkq6aY0dS7vIsO+u0ghonT2YkzFX1t6ra+Vi95wGMNLHfKPAh0B9dUBiXbmWrkPs6WVnJoJTLFY+2tTtzdcvnlqXepfrLaz9qqRL7whgzvwbAUyHs11lRfoiFD3+M8qTg+jpZ1aTRmHHuyaGLZZtewxMN+3v4ivAtqWvC20dPPpr5zlnneTf8mY3AYS4idSLyUoZ/lUnb3AqgDcAj3eyjWkTqRaQ+Fov1vXpHRPUhFgxyM30d38bL3k5239yLcUrRyeUb/1+Dtcm4Ht28J2Xu9RnnDvP+dv3eiKr2vlWQHYl8EcACAFeo6tHetq+oqND6+nojx3ZFac363jdyhG9BLiJbVbUihP1m1deAn73daWPjQcxbdfISxcoLh+P/5vnO6Cca9qdMonVqf0HjHTPyWkM+Be1tU1ezTANwC4Crgja8j6ISkGedVtT7RsS+zmBK+bCU8fPHt7+Z1/Hz5tiRLrMhLv+C8b/hkWRqzPynAM4A8LSINIjIckP7jRzXA/19pwi2fnea7TKign2dQdWk0fjJnPGJ5WWbXstboC98OPUVz01Tzy7IyxAzMXU1y9mqOkpVx8f/LTCx36hyOdC3f8/fl6Omsa+79y/jR+BDg09LLC/b9Fqo158famnF3J89j1f+2pJYt3DyGHx16tjQjhk1vAM0JK4FejHcq4mi7cFrJiJ5wO4zy/4UyuPmDrW0ovKnf8SzzW8n1t0567yCvZ68OwzzELkSngLgFUdqIX+UlQzCLxdehqKkCbnuqG00eslic+wIrlzyLPa+815i3U1Tzy74K1cyYZiHzHagTy4bgtcY5BSSCaMH45cLLkuZYfGG1Q244ke/R3PsSPdfGMATDftxxY+fwYHD/0is49BK9xjmeWAr0Hcvnun1A2zJDRNGD8aGmz6OgcUn46T5rfdwxY+fyfksfUldU5erVtYsvIxDKz1gmOfJ7sUzceEH8zdJvu1XBFRYykoGYd0NH0PZ0NNS1t+wugGX/KAu8JujzbEj+PTdf0i5IQjoGCOfMHqwsXp9ZOymoWz5fGNFb8K8uejOWecV5HhiWDcN5aKQexvouDsz/dFtADCouAgiwE+qJiQuJzzU0opb17yAp3b8FUX9gPZ2IDmRBp/WH/fPm1jQQR60t/vnoxhK1XnWbDLUx5acHtlpBcgvVZNGY+zwM7Hw4a14u+UY2to71h9p7ZimtvMO0lP6C9rbgePtHfF9Ir6dACjqJ7i98tyCPDHJFcPcIhOhvmbhZQV91kJumjB6MDbfOhWHWlrxjcca8PumGAYVFyUCHQCOtaWOChT1A07t3w9L/+0i3giUA4a5A5LHtw+1tGLC7U93u+1NU8/mu/kUGUMGFuOBL01MLG9sPIhrH9mKo8fbcUp/QZEIzh42CPfM+UjBPHg5LAxzxwwZWMw3L8lbU8qH4eXbp9suw0u8moWIyAMMcyIiDzDMiYg8wDAnIvIAw5yIyAMMcyIiDzDMiYg8wDAnIvIAw5yIyAMMcyIiDzDMiYg8wDAnIvIAw5yIyAMMcyIiDzDMiYg8wDAnIvIAw5yIyAMMcyIiDxgNcxH5hoioiAw1uV8im9jXFAXGwlxERgH4FIDXTe2TyDb2NUWFyTPzewB8E4Aa3CeRbexrigQjYS4iVwHYr6rbTeyPyAXsa4qS/kE3FJE6AMMzfOpWAN8C8OkA+6gGUB1fPCYiLwU9fp4MBfCW7SLSsKZgxubyRSb6Or4f9nb2WFMwgXpbVPv26lFEzgfwOwBH46tGAjgAYKKqvtnD19WrakWfDm4YawqmEGrKta/DqMUE1hRMlGsKfGbeHVV9EcCwpAPvBlChqq79dSMKjH1NUcPrzImIPNDnM/N0qloacNMVpo9tAGsKpuBqyqKvgQL8+eSINQUTqKY+j5kTEZF9HGYhIvKAE2Hu0u3SIvJDEWkUkRdEZK2IvN9iLdNEpElEdolIja06kuoZJSJ/EJGdIrJDRG60XVMnESkSkb+IyJO2a0nG3s5Yh1N9Dbjb29n0tfUwd/B26acBnKeqFwB4BcAiG0WISBGApQCmAxgHoEpExtmoJUkbgJtV9cMALgFwnQM1dboRwE7bRSRjb3flaF8D7vZ24L62HuZw7HZpVf2tqrbFF59Hx/XFNkwEsEtVX1XVVgCrAVRaqgUAoKpvqOq2+Md/R0eTjbBZEwCIyEgAMwGstF1LGvZ2V871NeBmb2fb11bDPAK3S18D4ClLxx4BYG/S8j44EJydRKQUwEcAbLZcCgDci47QbLdcRwJ7u1tO9zXgVG/fiyz62vilielM3S5tUk81qerj8W1uRcdLr0fyWVsSybDOiTM8ERkE4NcAvqaq71qu5UoAB1V1q4hMyfOx2dvZc7avAXd6O5e+Dj3MVXVqpvXx26XHANguIkDHS75tItLr7dJh1ZRU2xcBXAngCrV37eY+AKOSljtvJ7dKRAago9kfUdU1tusBcDmAq0RkBoBTAZwpIg+r6hfCPjB7OydO9jXgXG9n39eq6sQ/ALsBDHWgjmkAXgZQYrmO/gBeRUcoFAPYDuBcyzUJgIcA3Gv7/6mb+qYAeNJ2HRnqYm+frMG5vo7X5WxvB+1rF94Adc1PAZwB4GkRaRCR5TaK0I43qq4HsAEdb8Y8pqo7bNSS5HIAcwF8Mv6zaYifOVA0WO9tR/sa8KC3eQcoEZEHeGZOROQBhjkRkQcY5kREHmCYExF5gGFOROQBhjkRkQcY5kREHmCYExF5gGFOROQBhjkRkQcChXlvj3kSkSkicjhpToN/N18qkVnsa/JJr1PgJj3m6VPomL5yi4isU9WX0zZ9VlWvDKFGIuPY1+SbIGfmTj7miaiP2NfklSBhHvQxT5eKyHYReUpEzjVSHVF42NfklSBPGgrymKdtAEar6pH4HMD/DeCcLjsSqQZQDQADBw68qLy8PLtqI+zF/Ydx/oj3eXs827Zu3fqWqpZk8SXG+hoo7N7u9OL+wwBgvO8KrZfTBe3tIGHe62OeNOlZeapaKyL3ichQVX0rbbsVAFYAQEVFhdbX1wc4vB9Ka9ajfvFMb49nm4jsyfJLjPV1/PMF29udSmvWA4Dxviu0Xk4XtLeDDLNsAXCOiIwRkWIAcwCsSzvYcIk/7FBEJsb3+3Z2JRPlFfuavNLrmbmqtolI52OeigA8oKo7RGRB/PPLAcwGsFBE2gC8B2CO8hFG5DD2NfkmyDALVLUWQG3auuVJH/8UHc8XJIoM9jX5hHeAEhF5gGFOROQBhnkedb7b78txiMgdDHMiIg8wzInIKJOvDPkqMziGORGRBxjmREQeYJgTEXmAYU5EzuPYee8Y5kREHmCYExF5gGFOROQBhjkRkQcY5kREHmCYExF5gGFOROQBhjkRkQcY5kREHmCYExF5gGFOROQBhjkRkQcChbmITBORJhHZJSI1PWx3sYicEJHZ5kokCgf7mnzSa5iLSBGApQCmAxgHoEpExnWz3V0ANpguksg09jX5JsiZ+UQAu1T1VVVtBbAaQGWG7W4A8GsABw3W552wp/LkVKGBsa/JK0HCfASAvUnL++LrEkRkBIBZAJb3tCMRqRaRehGpj8Vi2dZKZJKxvo5vy94mq4KEuWRYp2nL9wK4RVVP9LQjVV2hqhWqWlFSUhKwRKJQGOtrgL1N9vUPsM0+AKOSlkcCOJC2TQWA1SICAEMBzBCRNlX9bxNFEoWAfU1eCRLmWwCcIyJjAOwHMAfAvyVvoKpjOj8WkVUAnmTDk+PY1+SVXodZVLUNwPXoeDd/J4DHVHWHiCwQkQVhF0gUBva1WelvvJt4Iz6MffosyJk5VLUWQG3auoxvCqnqvL6XRRQ+9jX5hHeAEhF5gGFuQVgvF/kylKhwMcyJiDzAMCci8gDDPA9sDX9w2IVsYv/lF8PcE/zFIZ+wn7PHMLeEzUpEJjHMiYg8wDC3yNTZOc/yiYhhTkRO4clJbhjmEcfGJyKAYW4dw5iijj3sBoa5A3L5ZSitWc9fIio47PnuMcwdkU04s6EpKtir+RNoClzKn+Tm3714Jn8ZiCgQnpk7jEFOREExzIkoVNmclPAEJncMcyIiDzDMiYg8wDAnotAFGT7hEEvfBApzEZkmIk0isktEajJ8vlJEXhCRBhGpF5GPmi+VyCz2Nfmk10sTRaQIwFIAnwKwD8AWEVmnqi8nbfY7AOtUVUXkAgCPASgPo2AiE9jX5JsgZ+YTAexS1VdVtRXAagCVyRuo6hFV1fjiQAAKAmD/paPt4zuMfZ1nPfUi+7TvgoT5CAB7k5b3xdelEJFZItIIYD2Aa8yURxQa9rUjGORmBAlzybCuyxmKqq5V1XIA/wrg9ow7EqmOjz3Wx2KxrAolMsxYXwPs7aDSp61gkJsT5Hb+fQBGJS2PBHCgu41VdZOIlInIUFV9K+1zKwCsAICKigq+ZCWbjPV1/PPs7SwwxM0Lcma+BcA5IjJGRIoBzAGwLnkDETlbRCT+8QQAxQDeNl0skUHsa/JKr2fmqtomItcD2ACgCMADqrpDRBbEP78cwP8GcLWIHAfwHoD/k/TGEZFz2Nfkm0CzJqpqLYDatHXLkz6+C8BdZksjChf7uu9sDJeU1qzH7sUz835c1/EOUCIiDzDMiYg8wDAnIvIAw5yIyAMM8xC5ci2tK3UQUXgY5kSUE54kuIVhTkSRwz8kXTHMiYg8wDAnIvIAw5yIyAMMcyLKGses3cMwJ6JI4h+UVAxzIiIPMMxD4tpZg2v1EJFZDHMiIg8wzImIPMAwJyLyAMOciMgDDHMiyopLb6a7VIttDHMiIg8wzEPg6tmCq3URUd8FCnMRmSYiTSKyS0RqMnz+8yLyQvzfn0TkQvOlEpnFvs6eiycELtZkQ69hLiJFAJYCmA5gHIAqERmXttlrAD6uqhcAuB3ACtOFEpnEvibfBDkznwhgl6q+qqqtAFYDqEzeQFX/pKrvxBefBzDSbJlExrGvyStBwnwEgL1Jy/vi67rzZQBPZfqEiFSLSL2I1MdiseBVEplnrK8B9jbZFyTMJcM6zbihyCfQ0fS3ZPq8qq5Q1QpVrSgpKQleJRnD8cUEY30NsLfJviBhvg/AqKTlkQAOpG8kIhcAWAmgUlXfNlNe9DAsI4N9nSWXe9vl2vIlSJhvAXCOiIwRkWIAcwCsS95ARD4EYA2Auar6ivkyiYxjX5NXeg1zVW0DcD2ADQB2AnhMVXeIyAIRWRDf7N8BnAXgPhFpEJH60ComMoB9nZ0onPlGocYw9Q+ykarWAqhNW7c86eP5AOabLY0oXOxr8gnvADUoKmcGUamT3BClfolSraYxzImIPMAwJyLyAMOciLoVxWGLKNZsAsO8QBVqw1NwUe6RKNeeK4a5IYXYPOQv9nP0MMyJyEuF9geJYU5UIIKGm08hWEjfM8PcgKg2QlTrptz19H9eWrPey57o7fvy5XsOdAcoEfnDl/DKlu/fN8/MiYg8wDDvo6j/tY96/UTUgWFOROQBhjkRkQcY5n3gyxCFL98HUSFjmBMReYBhTkTkAYZ5jnwbmvDt+yEqNAxzIiIPMMyJiDzAMM+Br0MSvn5fRIUgUJiLyDQRaRKRXSJSk+Hz5SLynIgcE5FvmC+TyDz2Nfmk1zAXkSIASwFMBzAOQJWIjEvb7BCArwL4kfEKHeP72avv318n9jX5JsiZ+UQAu1T1VVVtBbAaQGXyBqp6UFW3ADgeQo1EYWBfk1eChPkIAHuTlvfF12VNRKpFpF5E6mOxWC67sKpQzloL5Ps01tdA9Huboi9ImEuGdZrLwVR1hapWqGpFSUlJLrsgMsVYXwPsbbIvSJjvAzAqaXkkgAPhlOOuAjlbTSiA75d9TV4JEuZbAJwjImNEpBjAHADrwi2LXOB5oLOvySu9PjZOVdtE5HoAGwAUAXhAVXeIyIL455eLyHAA9QDOBNAuIl8DME5V3w2v9PzxPNQKEvuafBPoGaCqWgugNm3d8qSP30THy1TvFHqQl9asx+7FM22XEYpC7mvyD+8AJSLyAMO8B4V+Vt6JPwci9zHMu8EAS8WfB5HbGOYZMLgy48+FyF0McyIiDzDM0/Dss2f8+RC5iWGehEEVDH9ORO5hmMcxoLLDnxeRWxjmYDARUfQVdJiX1qxnkPcBf3ZE7ijYMGcQmcE/iERuKMgwZ/iYx58pkV0FF+YMnfDwZ0tkT6BZE33AoMmPzp+zrzMtErnK+zBniNvBUCfKL6+HWRjk9vH/gCg/vA5zIqJCwTAnIvIAw5yIyAMMcyIiDwQKcxGZJiJNIrJLRGoyfF5EZEn88y+IyATzpRKZxb4mn/Qa5iJSBGApgOkAxgGoEpFxaZtNB3BO/F81gGWG6yQyin1NvglyZj4RwC5VfVVVWwGsBlCZtk0lgIe0w/MA3i8i/8twrUQmsa/JK0HCfASAvUnL++Lrst2GyCXsa/JKkDtAJcM6zWEbiEg1Ol6uAsAxEXkpwPHzaSiAt2wXkSbyNcldIVZy0tgstzfW1wB7O0dO1RTvU6dqigvU20HCfB+AUUnLIwEcyGEbqOoKACsAQETqVbUiSJH5wpqCcbWmLL/EWF8D7O1csKZggvZ2kGGWLQDOEZExIlIMYA6AdWnbrANwdfzd/0sAHFbVN7KqmCi/2NfklV7PzFW1TUSuB7ABQBGAB1R1h4gsiH9+OYBaADMA7AJwFMCXwiuZqO/Y1+SbQLMmqmotOho7ed3ypI8VwHVZHntFltvnA2sKxouaQurrnGrJA9YUTGRrko5+JSKiKOPt/EREHnAizEXkGyKiIjLUgVp+KCKN8du314rI+y3W0uPt5hbqGSUifxCRnSKyQ0RutF1TJxEpEpG/iMiTtmtJxt7OWIdTfQ2429vZ9LX1MBeRUQA+BeB127XEPQ3gPFW9AMArABbZKCLg7eb51gbgZlX9MIBLAFznQE2dbgSw03YRydjbXTna14C7vR24r62HOYB7AHwT3dyMkW+q+ltVbYsvPo+Oa4ttCHK7eV6p6huqui3+8d/R0WTW74gUkZEAZgJYabuWNOztrpzra8DN3s62r62GuYhcBWC/qm63WUcPrgHwlKVjO30ruYiUAvgIgM2WSwGAe9ERmu2W60hgb3fL6b4GnOrte5FFX4f+QGcRqQMwPMOnbgXwLQCfDruGdD3VpKqPx7e5FR0vvR7JZ21JAt9Knm8iMgjArwF8TVXftVzLlQAOqupWEZmS52Ozt7PnbF8D7vR2Ln0depir6tRM60XkfABjAGwXEaDjJd82EZmoqm/aqCmpti8CuBLAFWrv2s3At5Lnk4gMQEezP6Kqa2zXA+ByAFeJyAwApwI4U0QeVtUvhH1g9nZOnOxrwLnezr6vVdWJfwB2AxjqQB3TALwMoMRyHf0BvIqOUCgGsB3AuZZrEgAPAbjX9v9TN/VNAfCk7Toy1MXePlmDc30dr8vZ3g7a1y68AeqanwI4A8DTItIgIst7+4IwaMcbVZ23m+8E8Jiq7rBRS5LLAcwF8Mn4z6YhfuZA0WC9tx3ta8CD3uYdoEREHuCZORGRBxjmREQeYJgTEXmAYU5E5AGGORGRBxjmREQeYJgTEXmAYU5E5IH/D/1Jp8+D/9kwAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(2, 2)\n",
    "axs[0, 0].scatter(x_list, v_list, s=1)\n",
    "axs[0, 0].set_xlim(-4, 4)\n",
    "axs[0, 0].set_ylim(-4, 4)\n",
    "axs[0, 0].set_aspect(1)\n",
    "\n",
    "choose = (vel_eta_list >= -0.01) & (vel_eta_list <= 0.01)\n",
    "axs[0, 1].scatter(x_list[choose], v_list[choose], s=1)\n",
    "axs[0, 1].set_xlim(-4, 4)\n",
    "axs[0, 1].set_ylim(-4, 4)\n",
    "axs[0, 1].set_aspect(1)\n",
    "\n",
    "axs[1, 0].hist(v_list, density=True, bins=100)\n",
    "axs[1, 0].set_xlim(-4, 4)\n",
    "axs[1, 0].set_ylim(0, 0.5)\n",
    "axs[1, 0].set_box_aspect(1)\n",
    "\n",
    "axs[1, 1].hist(x_list, density=True, bins=100)\n",
    "axs[1, 1].set_xlim(-4, 4)\n",
    "axs[1, 1].set_ylim(0, 0.5)\n",
    "axs[1, 1].set_box_aspect(1)\n",
    "fig.set_size_inches(6, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAokAAADUCAYAAADwWLkhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABADUlEQVR4nO3dd3xUVf7/8ddnJh0IoSM1KL0JCAgiIlVQEHFF0bVXvrKi7rpiXevq6rq6PxELsna/6voFFRdYJYiuNAFFYJFiRJoIEoqUhEySOb8/EmbupENm5kz5PH3ch3dm7tz7zpBPcnLPveeIMQallFJKKaWcXLYDKKWUUkqpyKONRKWUUkopVYY2EpVSSimlVBnaSFRKKaWUUmVoI1EppZRSSpWhjUSllFJKKVWG9UaiiLhFZJWI/Mt2FqWihYiMFJGNIpItIneV8/rZIvKriHxbsvzJRk6lopXWmFKQYDsAcCuwHki3HUSpaCAibmAaMBzYAawQkdnGmO9KbfqlMWZ02AMqFeW0xpQqZvVMooi0AM4DZtjMoVSU6QtkG2M2G2M8wLvAWMuZlIolWmNKYb+7+e/AnYDXcg6loklzYLvj8Y6S50rrLyKrRWSeiHQJTzSlYoLWmFJY7G4WkdHAL8aYr0Xk7Eq2uxG4EaBWrVqndejQMTwBlarEN998nWOMaWTp8FLOc6Xn1/wGaG2MOSwi5wIfAu3K3VmpGuvYUWvsmOxfDgc8btu4dkj2Hcz9xoqvv46NGtP6UpGouvVl85rEAcD5JcWVAqSLyFvGmMudGxljpgPTAU47rbdZ/NXK8CdVqpTURNlq8fA7gJaOxy2Anc4NjDEHHetzReR5EWlojMkpvTNnjfXu3dusXKk1dsyYqYsCHn98y5kh2Xcw9xsrRGKjxrS+VCSqbn1Z6242xtxtjGlhjMkEJgCflW4gKqXKtQJoJyJtRCSJ4vqZ7dxARJqKiJSs96W41veGPalS0UlrTCki4+5mpdRxMMYUisjvgE8AN/CKMWadiEwsef1F4CLgf0SkEMgDJhhjSneXKaXKoTWmVLGIaCQaYz4HPrccQ6moYYyZC8wt9dyLjvXngOfCnUupWKE1ppT9u5uVUkoppVQE0kaiUkoppZQqQxuJSimllFKqDG0kKqWUUkqpMrSRqJRSSimlytBGolJKKaWUKkMbiUoppZRSqgxtJCqllFJKqTK0kaiUUkoppcrQRqJSSimllCpDG4lKKaWUUqoMbSQqpZRSSqkytJGolFJKKaXK0EaiUkoppZQqQxuJSimllFKqDG0kKqWUUkqpMrSRqJRSSimlytBGolJKKaWUKkMbiUoppZRSqowE2wGUUkopFfvGTF1U5rmPbzkz5o8dzbSRGGO8xgQ89hR6fevOl0T860lu/wlll8vxglJKKaXilnY3K6WUUkqpMrSRqJRSSimlytDu5ijl9fr7jo/kF/nWi0p1Nye6/d3Hgn/d2S19qKDQv02p3uZayf5vEbd2RSullFJxQ88kKhWFRGSkiGwUkWwRuauS7fqISJGIXBTOfEpFO60xpbSRqFTUERE3MA0YBXQGLhWRzhVs9wTwSXgTKhXdtMaUKqaNRKWiT18g2xiz2RjjAd4Fxpaz3S3ATOCXcIaLZ4WFheTm5tqOoWpOayxC5eXlUVhYWPWGKii0kRhF8gu8vuXXvALf4nLhW5ITXAGLS8S/OLZzu8S3JCW4fEuCK3A5mFfgW/I8Rb5FWdUc2O54vKPkOR8RaQ6MA16samcicqOIrBSRlXv27Alq0Hjypz/9iZSUFGrVqkWvXr2YNWuW7UjqxAWtxrS+giM/P5+uXbuSlpZGeno6EyZMIDs723asmKeNRKWiT3l3EJlSj/8OTDHGVNmiN8ZMN8b0Nsb0btSoUTDyxaUzzjiD2267jQcffJD8/Hx+85vfcNttt2FM6X8aFQWCVmNaX8GRnJzM2LFjefTRR7nmmmuYM2cOp512GgsWLLAdLabp3c1KRZ8dQEvH4xbAzlLb9AbeleLb1RsC54pIoTHmw7AkjBNz5szhs88+46mnnmLkyJGMHDkSgPvuu48777yTTp06IaWHDFDRQGssAhhjuOmmm/jtb3/LoEGD+POf/+x7bcqUKdxyyy20bNmykj2omtJGYoRzdu0615MTHbOkVPJL6Hh/P5XePinBf5yjBf7jO4fQcQ6To8JiBdBORNoAPwETgMucGxhj2hxbF5HXgH/pL6/g2rZtG1dccQWZmZnk5+eTkpLie83tdvO3v/3N99jj8ZCUlGQjpjoxWmMR4KWXXuLll1+mTZs2DBo0KOC1Vq1a8dFHH/kea42FhnY3KxVljDGFwO8ovqNyPfBPY8w6EZkoIhPtposPxhiuu+46CgoKeO+99wIaiKXNnj2bjh07sn///jAmVDWhNWbf5s2buf322znnnHOYMmVKhdsZY7j22muZOFH/WULBWiNRRFqKyEIRWS8i60TkVltZlIo2xpi5xpj2xphTjDF/LnnuRWNMmYvojTFXG2P+L/wpY9eHH35IVlYWjz/+OO3atat021atWrF161YeeuihMKVTwaA1Ztcdd9yB2+3mH//4By5XxU0VEaFhw4a8+uqrrFixIowJ44PNM4mFwB+MMZ2AfsCk8sahikfOu4iP5PsX513I4vivMsZUvVTGeZxEt8u3HC3w+pYj+YW+RalY5/V6+eMf/0jXrl2rdfaiR48e3HDDDTz33HNs3LgxDAmVim6LFy/mgw8+4N5776V58+ZVbn/ffffRpEkTvVEsBKw1Eo0xPxtjvilZP0TxKf2qvxuUUsoil8vFxx9/zOuvv05CQvWux33kkUdITk7mscceC3E6paJfv379mDVrFrfffnu1tk9PT+fBBx9kyZIlerdzkEXENYkikgn0BL4q5zX/GFM5OsaUUsq+Tp060atXr2pv36hRI2666Sbefvttdu4sfZOsUsrJ7XYzbty4Sq/1Le2aa66hWbNmPPXUUyFMFn+sNxJFpDbFI9bfZow5WPr1gDGmGsbmGFMFRd6A5WBeoW9xDowtgm8xzv+q2XVcHZV1RTuPn+R2+RZnl7in0BuwKBVLPvzwQy666CL27t173O+98847+fLLL2nWrFkIkikVG6677jqeeeaZ435fcnIy//znP3njjTdCkCp+WW0kikgixQ3Et40xOj2BUiqiPf3006xatYqMjIzjfm/Tpk3p379/8EMpFSN+/PFHXn31Vfbt23dC7x8wYACNGzcOcqr4ZvPuZgH+Aaw3xjxtK4dSSlXHoV1b+PLLL5k4cSJut/uE9pGXl8f111/P66+/HuR0SkW/GTNmICLcdNNNJ7yP5cuXM2zYsBNuaKpANs8kDgCuAIaIyLcly7kW81iTc8gTsCQnunwLgm+pqBvY2fVc5K3ZUvo/r/EvAcd35HLedV36awlmV7hSNm1b8jGJiYlcddVVJ7yPlJQUVq5cyTPPPKN3YSrlUFBQwCuvvMK5555LixYtTng/KSkpLFiwgNdeey144eKYzbubFxljxBjT3RjTo2SZayuPUkpVpKggnx0rPuGCCy6oUXeWiDBx4kRWr17NgW3rg5hQqeg2Z84cdu3axY033lij/XTv3p0zzjiD6dOn6x9iQWD9xhWllIp03gIPmWdewOTJk2u8rwkTJpCUlMRPK+cHIZlSsaF58+bcfPPNjBo1qsb7uvrqq9m4cSPffPNNEJLFN510VymlqpCYVocO513PmWeeWeN9ZWRkMGbMGD7+dAGdL5iEy60/hlX8GjN1kf9Bx0sZ98IyPr7lxOtszNRFeHJb4HInMv6OJ+hyYc3/sItn+tPJkpxDHt96gitw1pSAR46z5cbxgqlgdBlDzU6vi6lkBhfnS47DOOMnJQSenN57ON+33rBOco2yKWVDQd4R9v3wLY069g34hXa8v8ic7/2lcX+a9YQiz1FcqbWDllWpaHRg2wZcCUmkNzs5KPtLSqtDm8EXU6tRy6DsL55pd7NSSlVi99ovWTH9Ln7dHrwp9Rp3Op2uF91GojYQlWLDxy+xcsY9Qb2GsNP5E2nV/7yg7S9eaSNRKaUq8dM3C0it35SMzC5B3a/xFrHvhzUUFXiq3lipGJV/cB85m76h+WnDKB4ZL3gK8g6zf+t3Qd1nvNFGYhgVFHp9S25+oW9JcEvA4hQwHI3XvxQZx1LJkDaFRVUvlQ6L4ziO8/gVzfjidgUuuflFviW/wOtblIoGnsMHyNmwgma9hgb9F9iejStZ8v8mkbPp66DuV6lo8vO3n4Px0qzX0KDve+0/n2bFS1Mw3qKg7zteaCNRKaUqsPu/izHeIk7qMTjo+27QticJyWnsWv1F0PetVLT4efUX1G6aSZ2T2gR93027nYnn8AH2bV4b9H3HC71xRSmlKrA3+1tSMhpTt2X7oO/bnZhE4y792f3fxRQVFZ3wLC5KRasiz1EObFtPm7N+E/B8wB3PJcq7Uay87Zwad+6HKyGJXWv+Q4O2PWqUNV5pIzGMfv71qG89JaniXwjOS3eddzE771z2OjbyOh4UlZ2Opdz9BnScOR64S3WpuRy3LjvvYg64C7qSXjjn17n7oP/rb9UgreI3KWVJ6V86p152N0d/3RP0ruZjmp46iJ3fLGDg718I+CVWkyFAlIoW7qQUhj/6Id7Cgiq3rapBWJ6ElDQadujNrtX/ofO4W0JWx7FMu5uVUqoC4nKRWq9JyPbfuPPpuBKS2P3fxSE7hlKRLCE5jaRadUO2/5NOPYu8/bs5tPOHkB0jlumZRKWUKsfGOTMoKvTQeezNITtGQnIaA25/njonBWd8OKWihbewgK9euIOTh1xCky5nhOw4TU8dREZmF+o0zQzZMWKZNhJDLM/jv6vKU+jvO05zdMOWHRqq/G5l5xhShc4u5iL/ekFR4M4Kvf5jOlYDuo7djjuqk9yBJ5edrzkH/ZaArufS+cs/TkGhf8PcfP/nkpYcv9diichM4BVgnjEVDZGuws14vWxb+i/qn9wt5Meq27JDyI8Rz7TGItPeH1az9/tvaHP2+JAeJzG1to5HWgPa3ayUXS8AlwHfi8hfRKSj7UAKft2+kfyDe2nSdUDIj+UtKmT97Bd1LufQ0RqLQLvXLsKVmESjDr1DfqyDOzez6o1HyD+4L+THijXaSFTKImNMljHmt0AvYAswX0SWiMg1IpJoN1382r1uCYiLxl36h/xYLncCu/+7mO1fzQ35seKR1ljkMcaw+79LaNi+N+6klNAfr6iAn1Z+yp4Ny0N+rFij3c0htuuA/47etCTnx+3vhy3dW+vsFq6oi7mgyL+RxzE49eH8woB97cn1z5180OO/g6y2I0uTNH+R1koO/JZITvT/HWEcXdHOKZpdrurdMZbqvNPZcad3m8a1qvX+WCUiDYDLgSuAVcDbwJnAVcDZ9pLFr5wNK8lo1TGkF9Q7Ne7cjx+/+D8K83NJSNY7/4NNayyy5Ob8RN6+nzll6KVhOV5683Ykp9fnl++W0qLvyLAcM1bomUSlLBKRWcCXQBowxhhzvjHmPWPMLUCFF9KIyEgR2Sgi2SJyVzmvjxWRNSLyrYisFBEdU6WajDHUaX5KSGaAqEjjzv0xRYXkbNTZV4JNayzyFHmO0rhzPxqGoasZikcpaNypH3vWL8dbVFj1G5SPnklUyq4ZxpiAfkYRSTbG5Btjyv0JKiJuYBowHNgBrBCR2cYY5ySlC4DZxhgjIt2BfwJ6LVY1iAjdL7mjyu1Kj9tWemzD4xnXrf7J3UhITuOX75bRtPvAar9PVYvWWIRJb96WvhP/GtZjNurcj+1fzeXAlu+of0r3sB47mumZRKXserSc55ZW8Z6+QLYxZrMxxgO8C4x1bmCMOWz81yrUouxVDaoC+Qf3BVzmEQ6uhERO6nE2rgT9uz0EtMYiiPEWkX9of9iP26hjH2o3aU1B3qGwHzua6U+kEHDOgHLoqP/UdmrAUC/+bYwJvKbP6/gFVVTBdYjOoXX25/qvNVy+M7D4Hnj6E/+DHY4/gk9q51u9/w+jfev9W2QEvL9+WpL/gWNVHH9fOH8yVnZ5osvxJ8lhx+fi/Lyqe31jtBORpkBzIFVEeuK/SDWd4m6xyjQHtjse7wBOL+cY44DHgcbAeTXNHC+WTp1MnWancNo1D4X1uKf+9u6wHi/WaY1Fpl93ZLPoqevpc8NfaNIt9KMHHJOYWpuz730rbMeLFdpIVMqOc4CrgRbA047nDwH3VPHe8lrSZUfbNOYD4AMROQt4BBhW7s5EbgRuBGjVqlVVuWNa3oE9HN69lZb97P2+LyrIr3ojVR0RUWNaX4FyNq4EIKN1JyvHN94ijDG43Nr8qQ79lJSywBjzOvC6iPzGGDPzON++A2jpeNwC2FnJsf4jIqeISENjTE45r08HpgP07t07rrvM9m4qvnHkRC6oP5G5ZUtb/uIfi0eq/334bpqJVZFSY1pfgXI2raROs5NJTq8f9mMf/mU7i5++iW4X30GzXkPCfvxopI3EEPjZMbxL7ZSqP2JvqeufnN2vhY4ZVPIdQ938mufvYl6yw9/F/NhLXwbsa/HUy33rnVuk+9Y37vRflzHgDv/Pz7tvCrz4fmBLfyG7KpgcPcnR9SyluosDZmZx/HFey/G5OD+v5vVSyz1GrBGRy40xbwGZIvL70q8bY54u523HrADaiUgb4CdgAsWDBTv33xb4oeSi+l4UXyywN2hfQIzK2fQ1SbXqkt7sFCvHT63flB3LP6GgoIDERB3Crya0xiJPkSeffT+sofWZF1g5flqDkzDGkLNppTYSq0kbiUrZcWxwyOOeL8oYUygivwM+AdzAK8aYdSIyseT1F4HfAFeKSAGQB1xiwn03RpQxxrBn40oadjgNcdm5p69hhz5sXfQhy5YtY+BAvcu5hrTGIsz+H9fiLfSEbeib0lzuBBq068meDSvDfnNatNJGolIWGGNeKvn/Cd0dUTKkx9xSz73oWH8CeKImGeONMYbuF99BYq30qjcOkYbteoK4yMrK0kZiDWmNRZ46J51Mt0v+SIO2p1rL0Kh9b3av+ZLcnAqvHlAO2kgMgUN5/jt3nTOYOLtrK/sjpsg5y0oFdzRn/3rYt/7Ua1/51ldPuyRgXxV133ZoVse3vvb5Cb71026fFbDdSZP93c9dE/yzTyQ4upWLXI67k0tf7+34Op1dz4mO2Vucnxf1yo0bc0Tk2cpeN8ZMDlcWVczlcoX1bsvyJKbVIaNVR+bPn89DD4X37upYozUWeZLT69N6wPlWMzTsWHwWM2fjCuCSyjdW2khUyhKdWiPCvPfee/y6I5e6LdpVvXEInTxkAjcPaWM1Q4zQGosgBw4cYNvSOTTtNoCk2hnWctRq1JL2o64lo3VnaxmiiTYSlbKg5M5LFSEKCwu54YYbqNd9SLVmWwmlZj0Hc/nlOsNbTWmNRZbPPvuMNe/8hdpNplHfYiNRRGg/6hprx4822kgMgaOObuE6qf6P2NkR6+x69XoJ4OyKLnTc6Xwg3+Nb/3DVbt/6qBH+v4hO5O7gkzJSfOvnntMl4DXncVoM9O87JdE/MHiS47uodDd6wPX/jtecN0EfLSgi3ojI340xt4nIx5Q//prdPpk4s2LFCg4dOkS79qfZjgLA999/z65du/S6xBrQGossWVlZuJNTI+IMnreokP2b17J1a0tat25tO05E00aiUna8WfL/p6ymUEDxLzARoWH7XrajADB58mS2bNnC+vXrbUeJZlpjEWT+/Pk0aNszIgaxLsg9xNKpk3m76SHuuaeqcdXjm87drJQFxpivS/7/BcXzyO4H9gFLS55TYZSVlUWvXr1IqlW36o3DYPjw4WzYsIEdO3bYjhK1tMYix5YtW8jOzrY29E1pyXXqkd68LVlZWbajRDz7TfoY5HH0H0u5szsFqmCMaiBw7ub9R/3dzfNf8F9us3tppTfxHZfnL+oW8LhJf/8NgL/t459btnEtfxd1dUebChhY2/HAU+QtZ+v4ICLnAS8CP1B8RUIbEbnJGDPPbrL44fF4+Pbbb5k4cSLfVb15WAwbVjy72/z587nmGr1+qia0xuxbvnw5AI0ipJEIxbMqLV40i9zcXNLSqprKO37pmUSl7PobMNgYc7YxZhAwGHjGcqa4kpSUxK5du5gyZYrtKD7dunWjcePGLFiwwHaUWKA1ZtnFF1/Mzp07qd0003YUn4YdeuPxeFi0qObTacYyq41EERkpIhtFJFtE7rKZRSlLfjHGZDsebwZ+sRUmXqWmplK/fvjnkq2IiDBkyBAWLlyoM0PUnNZYBDjppJMCepBsq39ydxITE/nss89sR4lo1rqbRcQNTAOGUzyZ+goRmW2MiZQeH6VCRkQuLFldJyJzgX9S3HM/nuJ5Y1WYXHHFFQwdOpSrr77adpQAjz/+OHXq1ImoX6zRRGssMqxbt4577rmHJ56IrMlpEpJTWblyJZ06dbIdJaJVq5EoIinAzcCZFBfZIuAFY8zRGhy7L5BtjNlccox3gbEQMZcFnTDnjCmB1+E5t3LOvlLxmQLnS3mFjqFixH8SOCkheCeEy+zLcZyjQbx20PlZOIf5iSNjHOu7gUEl63uIm3ln7Nu9ezdvvfUWnTvbH5ajtMzMTNsRop3WWAT45JNPmD17NtOmTQNybMcJ0L17d9sRIl51zyS+ARwCppY8vpTi4QXG1+DYzYHtjsc7gNNLbyQiNwI3ArRs1aoGh1Mqchhj9G6ECHCsq+nYjSKR5tVXX2Xbtm088MADtqNEHa2xyJCVlUXHjh1p0aIFsMV2nAC//vorDz/8MOeeey5Dhw61HSciVbeR2MEY45yRe6GIrK7hscvrQylvwNPpwHSA007rHZennFTsKjlLfx3QBfDdMm6MudZaqDiSlZVFRkYGvXpFxviIpS1btox33nmHe++9l4QEHYziRGiN2ePxePjiiy+49trI/KjT0tJ4+eWXyc3N1UZiBarbT7lKRPodeyAipwOLa3jsHUBLx+MWwM4a7jMiuEV8izH4lkDGsVRMxL+kJrh9C8brWzyF/qWmnPvyFHoDjpPidvmWmnJ+Li4R3xKH3gSaAucAX1BcB4esJooTxhjmz5/PkCFDcLvdVb/BgmHDhnHo0CFWrNBL6GpAa8ySZcuWkZubG7Fn6hMTExk0aJCOl1iJ6v62Px1YIiJbRGQLxQOTDhKRtSKy5gSPvQJoJyJtRCQJmADMPsF9KRWt2hpj7geOlMw1ex7QrYr3qCA4cuQIPXr0YMyYMVVvbMngwYMBdCicmtEasyQ/P5++ffsyaNCgqje2ZNiwYWRnZ7N161bbUSJSdfsvRgb7wMaYQhH5HfAJ4AZeMcasC/ZxlIpwBSX/PyAiXYFdQKa9OPGjdu3azJ4d2X+XNmzYkJ49e7JgwQLuu+8+23GildaYJcOHD2f48OG2Y1TqWDfzggULIrZb3KZqNRKNMSFpYhtj5gJzQ7Fvm5Ic3bHG0Z1sjL87tbo9q84u2HopSb714f9zlW990sy1vvWXL3FeOnr8/uf9wBPDzuM4j+92Ve8LqOjGbefnkhLEu7Oj0HQRqQfcT/GZ9Nol6yrEDh48SHp6uu0YVRo1ahTLly/HGKPD4ZwYrTELPB4PXq+XlJSUqje2qEuXLnTt2pXDhw/bjhKR9EpopSwyxswoWf0CONlmlnhSVFREmzZtqN97NB1H32A7ThljpjpmgWg6ivnz/2wvTJTTGrPjk08+4eKLL2bp0qX06NHDdpwKiQhr166tesM4pY1EpSwSkQbAg8AAiu9i+hJ4xBiz12auWLdq1Sr27dtH65Pa2I5SbUVFRRF7g00k0xqzY/Lf3sBTZLjns324v4zMqe8C/hgDjLeIf90auddP2qCNxBBISfT/IPc6bjg2zp/vjm7Y0j1IzseJbv+DjGR/d+8FPZv41v/w/770re86p0PAvppmVH2q/6f9eb71eZ8GjmX+1OQzyz1+gqv8rvPKesOcPc/Oz8X5ecWhd4H/AL8pefxb4D0gMm8HjBHH7mZs2C4yh74p7dZbb+Wrr75i2bJltqNEI60xC3I2rqTBKafiTkyqemPLCvKO8OVT15N55jjQRmKAuL4YTKkIUN8Y84gx5seS5VEgw3aoWJeVlUX37t1JTo+c+Zor07BhQ5YvX05OTmTNWBEltMbCbOfOnRzetYWGHXrbjlItiam1ECBn09e2o0QcbSQqZddCEZkgIq6S5WJgju1QsSwvL49FixZF1eC5Q4cOxRjDwoULbUeJRlpjYXZsyKZoaSQCNGx/GnuzV1FQUFD1xnFEu5tDID010bd+tMA/37JzXuSALtpS73c7XnQ7uptTk/zdsm3r1vat33GNfzbDLv/zTsC+lvztIt96h2Z1fOvf7TjoWz/rzlm+9XtuGhjw/nZ1/e9xdgs7c7kr6WN2vuTsYi5wzANdJzX+vg1F5BDFPfAC/B54q+QlF3AY0HnYQsQYw7Rp0+jZsycPLM61Hada+vTpQ506dViwYAHjx9dkNtT4oTVmz+mnn07HMRNJb3aK7SjV1rBDb7Yu/ogVK1Zwxhln2I4TMeLvt7NSEcAYU6fqrVQopKWlcd111xU/WByZF9SXpjNDHD+tMXvat29P2+G/tR3juDRo1wtEyMrK0kaig3Y3K2WZiJwvIk+VLKOr+Z6RIrJRRLJF5K5yXv+tiKwpWZaISM0G0IwhM2fOZPv27bZjHLfrr7+eiRMnUlRUVPXGKoDWWPj89NNPfPzxxxR5jtqOclySaqXT7pyr6N+/v+0oEUUbiUpZJCJ/AW4FvitZbi15rrL3uIFpwCigM3CpiHQutdmPwCBjTHfgEWB6sLNHo7179zJ+/HheffVV21GO29ixY7njjjt0GJzjpDUWXjNnzuT8888n/9A+21GOW4dzr4v4GWLCTbubQ6Bp3WTf+jrHtX9pyf4f7pXNnOByDC+T6GzHO8aQqeu47nFAC/8dmvf/z1kB++p30wz/g10/+Ndb+H/ePfwH/6yLfU7KCHh/uuN6wZREf5ZEx6wyzryuMl+XY8YZx/qRo4W+9cyGkT/rRQidC/QwxngBROR1YBVQ5syFQ18g2xizueQ97wJjKf4FCIAxZolj+2VAiyDnjkoLFizAGMOIESNsRzkhe/fuJTs7m9NPP73qjdUxWmNh9Omnn9K2bVvSGjSzHeW4GWP47rvvSEtLIzMz03aciKBnEpWyL8OxXrca2zcHnP2lO0qeq8h1wLzjjxV75s+fT926dendO3ruunSaPHkyY8eOxVQ036WqSIZjXWssRDweD59//nnU/hHmLfDQs2dPnnvuOdtRIoY2EpWy6zFglYi8VnKG4+uS5ypT3mnoclsNIjKY4l9gUyrcmciNIrJSRFbu2bOnmrGjjzGGTz/9lKFDh5KQEJ2dKMOGDWP37t2sW7fOdpRoYrXG4qW+AJYuXcqRI0eitsvWnZTMgAEDfEP4KO1uDgln92vtFP9H7BwChgpmLIFSLXfnA8cQOs4fYc4u3rNbNwjYV6cnL/WtH/b4u3jTk/zd1Y3S/N3jtZIDvyWSK+hiDpxxpeKvxRj/E0WOD6CW43Nxfl7xRERcgBfoB/Sh+F91ijFmVxVv3QG0dDxuAewsZ//dgRnAqMqmIDPGTKfkeqrevXvH7CmqrVu3sm3bNu6++27bUU7YsbEds7Ky6Nq1q+U0kS8Saixe6gvgP//5D263m8GDB/OPN6JzPuShQ4dy3333sWfPHho1amQ7jnV6JlEpS0qukfqdMeZnY8xsY8xH1fjlBbACaCcibUQkCZgAzHZuICKtgFnAFcaYTUEPH4UyMzPZtWsXl156adUbR6hWrVrRrl07HQqnmrTGwuvee+9l/fr11K1bnR79yDRsWPFsjTpwfTE9k6iUXfNF5A6K55I9cuxJY0yFtwYaYwpF5HfAJ4AbeMUYs05EJpa8/iLwJ6AB8HzJmd5CY0x0XogXRE2aNKl6owg0Zqp/PEdP48588cUCCgoKSExMrORdqoTWWJi4XC7atWtnO0aNnHbaaaSnp5OVlcXFF19sO4512kgMsSZ1U3zrPx/wjxuVlODooi11+YvLcX7XZZx3DvufL+5FKeac8cTZJQxQ29F97HV0dDj3leB2dikHZnHOrFJRF7NzX6Wvp3fe0ZyX7x/frXn9VBQA11J8rdPNpZ4/ubI3GWPmAnNLPfeiY/164PogZYx6BQUFTJgwgUmTJvHMuiTbcWrk5CGX8ukrT0btdZUWaI2FwSeffML777/Pk08+Sf360TEnenkSEhL497//TceOHW1HiQja3ayUXZ0pHo9tNfAtMBXoYjNQLPrqq6+YNWsWBw4csB2lxmo1ak779u0rHUZLBdAaC4OZM2fy/vvvk54e/UOa9e/fn3r16tmOERG0kaiUXa8DnYBnKf7l1ankORVEn376KS6XiyFDhtiOEhRz587l0UcftR0jWmiNhdixkQOGDBkSE2e4PR4Pf/nLX5g3L+5HNdLu5lBzDqCd6Ohidnb9luohDuh+loAuZue6/4Fb/Dtzu0v39/p37nzFeQ6ioq5jKD1QdvkZA8MHHt95R3dyov+zSE3SWSNKdDDGOKfzWigiq62liVHz58+nb9++ZGRk2I4SFF988QXPPPMMt99+O7Vq1bIdJ9JpjYVYdnY2W7du5c4777QdJSgSExN59tlnOeussxg1apTtOFbpmUSl7FolIv2OPRCR04HFFvPEnP3797N8+fKoHeC3PMOHD6egoEDvwKwerbEQmz9/PkDM1JiIMHz4cObPnx/3c6VrI1Epu04HlojIFhHZAiwFBonIWhFZYzdabNi9ezd9+vThnHPOsR0laAYOHEitWrWYO3du1RsrrbEQc7lcDBw4kFNOOcV2lKA599xz2bdvH8uXL7cdxSrtbg6jJun+O513H/Tf6ZzoDvxnqOh69Iq6oZ13NyeUmhSgOrN3SWXdyFLuaoVKHy/P4/8rrFlGCqqMkVVvomqiY8eOLFu2zHaMoEpOTmbYsGHMmTMHY4zexFI5rbEQmzhxIhMnTrQdI6hGjBiB2+1m7ty59O/f33Yca7SRqJRFxpittjPEMq/Xy3lPZ5GQnGY7SlCNmbqIrant2Zf/LTk5OTozRCW0xkLr8OHDpKWl4XLFVsdkvXr1OPvss9m3r8LhNONCbP2rKqWUw9dff82nd41mz8aVtqMEXat+ozn73re0gaisuueee2jXrl1MXrv36aefMm3aNNsxrNJGolIqZs2ZMwevt5C6zdvajhJ0UnLmxhswKbxS4WOMYc6cOXTq1Am3O/ZGrHBpjWl3czglJ/rb5GmOIWCKvIEX8rnczmFnKrgosILxbE5k9vgKh7OpJq/jQsTSX0stx4wviQn6N4kKrzlz5lCvdReSamfYjhISP32dRZMm49i0aVOFg/86p/QD+PiWM8MRTcWBjRs3snnzZu644w7bUULmkksuwev18v7779uOYoU2EpVSUcfZ8Cnd6Dn22tGDe1m5ciUdzrshrNnCKbVeY3Jychj2+6k061U8ULg2AlW4zJkzB4CZvzRmbqk/RqJV6T+q1vyUz77VC/B4PCQlRfeUnidCG4lKqahW+of6MXu+K76juXGX2L0zsV5mFxLT0vnlu6W+RmJFn4dSwTZnzhy6du1Kav0mtqOETJPO/dm2+CO+/PJLhg4dajtO2Gkj0ZIGtZN96z8fOBrwWoLz0o6KhqCpYCYWYwK7jgNfK//56nK+31TwvKcwsLu5Wb1klLKh/ik96HzBJNJj8HrEY8TlplGnvvzy3TKM1+u7TlGpcJgyZQpHjx5lxjbbSUKnQfteJCcnM3fu3LhsJOpPFKVUTKrVqDknD5kQ82MINulyBp7DBziw9TvbUVScOeeccxg7dqztGCGVkJzK4MGDmT17NqY6Aw/HGD2TqJSKOb9u30Te/l007nIGLnds/5hr3LkfJw+ZQFLt8m9cUSoU3n33XTp37kz37t1tRwm5SZMmsXnzZoqKikhIiO2fJ6XF11cbQZwnNxrWCbwYdv+RAt96WrK/77nSmVHK2abM4wq6nitjKrpf2vF0foF/eIBG6fF3Ya+KPFu+nMXP337OiMc+th0l5BLT6tD5gkm2Y6g4cvToUW644QYuu+wyXnrpJdtxQm706NG2I1hjpbtZRP4qIhtEZI2IfCAiGTZyKKVij7eokN1rFxWfRUxItB0nLLxFheRs+prcvT/bjqLiQFZWFocPH2bcuHG2o4TNvn37mDlzpu0YYWfrmsT5QFdjTHdgE3C3pRxKqRizb/NaPEd+5aRTz7IdJWwKcg+xbNrv2b5sju0oKg7MmjWL9PR0hgwZYjtK2Lz55ptcdNFFfP/997ajhJWVRqIx5lNjTGHJw2VACxs5IkVSgitgqZOS4Fs8hV7fYgy+RaR6S0Wq+36XiG9xyi/0+pb01ATfkuh2BSxKhduu1V/gSkyiUafTbUcJm+Q69WhwyqnsWvMf21FUjCssLGT27NmMGTMmrsYNPHbWdNasWZaThFck/Ba/FphX0YsicqOIrBSRlXty9oQxllIqGh3c8T2NOvYlITnVdpSwanrqWRz6+UcO747h8UiUdRs2bCAvLy+uupoBWrVqRZ8+fbSRGCwikiUi/y1nGevY5l6gEHi7ov0YY6YbY3obY3o3aqgT2SulKtf/1ufocfm9tmOEXdPuxd3rejZRhVLXrl3JyclhzJgxtqOE3YUXXsjy5cvZvn277ShhE7JGojFmmDGmaznLRwAichUwGviticfBhyqRluz2LckJLt9SWGR8i7PruTI16ZKG4nmZjy0FRV7fkpbk9i2pjkUpm4wxiAiJqbVtRwm71HqNyWjdmZxNX9uOomLUsV/VqampcdXVfMyFF14IwMKFCy0nCR9bdzePBKYA5xtjcm1kUErFFm9RIV88dgVbl8y2HcWa0659mL4T/2o7hopR8+bN49RTT+WHH36wHcWK9u3bs3nzZq688krbUcLG1jiJzwHJwPyS2RCWGWMmWsqilIoBOZu+5vDurSTVyrAdxZrUerE7h66y7+2332bHjh20bNnSdhRr2rRpYztCWNm6u7mtMaalMaZHyaINRKWOg4iMFJGNIpItIneV83pHEVkqIvkicoeNjOG28+ssElJr07hzP9tRrNr+1TyW/H0SxuutemNVIa2xQEeOHOHDDz9k/PjxcdnVfIzX6+Wyyy7j0UcftR0lLHTGlQhXO8X/T5TnKfKtO2c5SUrwt/VLX2PovGbR+VpF1zJ6S73gKfQfx3nNYUqiXn9oi4i4gWnAcGAHsEJEZhtjnJP37gMmAxeEP2H4FXmO8vPqL2jWcwjuxPj9BQYgLjf7Nq9h349raXDKqbbjRCWtsbJmz55Nbm4ua5K7MmbqIttxwqr01/vVqh9YunQp99xzDy5XJAwSEzqx/dUpFZv6AtnGmM3GGA/wLjDWuYEx5hdjzAqgoLwdxJrd/11MUX4ezXuPsB3FuqbdB+JOSuWnFZ/ajhLNtMZKefvtt0mp15j6J8f+XM1Vad5nBFu2bGHJkiW2o4ScNhKVij7NAecYDDtKnjshAWOR7onOsUhrN23DyYMvoUFbPXOWkJxK0+4D+XnVZxQV5NuOE62CVmOxUF8Al19+Oe1HXovE+Jmz6mjafSBpaWm8+eabtqOEnP5rRxHnUDN1UhN8S6HX61ucM7R4Cr0BQ9h4vfgXx/PO2VMCtvFC3dRE35KS6PYtyqryBi464WGkAsYibRSdY5GmNzuZzuN+h7j0exOg5emjKMg7zM/ffm47SrQKWo3FQn0BTJgwgVb9z7MdIyIkJKcxfvx43nnnHQ4fPmw7TkjpNYlKRZ8dgPP2whbATktZrNu15kuS6zagXuvOtqNEjAbtetF64IXUahy/d6HWkNZYCa/Xy7Rp07jkkktsR4kokyZNonnz5hQUxPbVBtpIVCr6rADaiUgb4CdgAnCZ3Uh2eIsKWfv+09Rt3lbHB3QQl4tu42+3HSOaaY2VyMrKYvLkyRSfBW1hO07E6NOnD3369LEdI+S0kRil3C5/b0h6aqJv3estdXdykf/uZOeNy27xX2mQmuTfl6uqKViUdcaYQhH5HfAJ4AZeMcasE5GJJa+/KCJNgZVAOuAVkduAzsaYg7Zyh8Ke774i/9ccWmmDqFyHd2/j8O6tNO0+0HaUqKI15vfyyy9Tv359xo0bx9vTV9iOE1GKioqYN28ebdu2pWPHjrbjhIQ2EpWKQsaYucDcUs+96FjfRRz82f/jF++TUrcRjbucYTtKRNo07x/sWb+chh1mkZCcajtOVNEag23btvHBBx9w++23k5ycbDtOxDl06BAXX3wxl112GTNmzLAdJyT0xhWlVFQ6+FM2OZu+JvOs3+By69+75ck86zcU5B1m+1dzq95YqVKmTp0KwC233GI5SWTKyMjgyiuv5K233mL37t2244SE/mSNMS5XYHdxit7tqWLU4d3bSK7bkFYDzrcdJWLVa9ONjMwu/Ljwn2SeeYHe/a2qbczURXyzcDVNTj2bSR9tA7bZjhRRjg2wfbjhWeTnv8TAq+9i07xXLacKPj2TqJSKSs16DWHog++TlFbHdpSIJSKcPPgScvfuZNfa+JolQ9Vcr6v+RI8r7rMdI6LVbtKKJl0HsPXLD8jNzbUdJ+i0kaiUijqHd2/FGKPdzNVw0qlnUbtpJrk5cTmCizoBHo+HI3t2AGiNVcPJQybgTk7lhx9+sB0l6LSRqJSKKjk5OXz51xvY8PFLtqNEBXG5OWvKq5wy9FLbUVSUeP311/n8z5dz8Kds21GiQv1TTmXw/e/QrVs321GCThuJSqmo8tRTT1FUcJQWfUfajhI1jp0NOrBtA8ac8OQ8Kg54PB4effRR6rbsQJ1mp9iOExVEBJc7gaNHj7JmzRrbcYJKG4lKqaixZ88ennvuOZr1Gkadppm240SV3WsXs+ipG5g7V+90VhV77bXX2LZtG+1HXYvouLnH5corr2TkyJHk5eXZjhI02khUSkWNhx9+mLy8PNqPvMp2lKjTqPPppDVoxv3334/X6636DSruHD58mAcffJB+/frRqFNf23Gizs0338zPP//M888/bztK0GgjUSkVFY4ePconn3zCxIkTqd2kte04UcflTqDD6BtYtWoVr7zyiu04KoKMmbqIMVMXMezO6ezesxd3/6v0LOIJ+NvaBBp1Op277nuAEX+e7RsmJ5ppI1EpFRVSUlJYs2YNjz/+uO0oUatZr6EMHDiQu+++m/3799uOoyJMo059GfrQ/1GvTVfbUaJWlwsnU1SQHzM31mkjUSkV8VatWkVubi4pKSmkp6fbjhO1RIRnn30WEWHdunW246gIYYxh3+biGy6S69SznCa61W7SijZnj+fwL9soKvDYjlNjOgCSUiqi7dmzh3POOYdBgwbx/vvv244T9Xr06MG2bdtISUmxHUVFiJ9WzufbNx+hzw1/oUm3AbbjRL0O516Py52AuKL/PJw2EpVSEe2WW27hwIED/OlPf7IdJWakpKTg9XqZMWMGl156KXXq6Kw18aL0dXJHD+5l3cy/Uy+zK4279LOUKra4E5MAyD+0n55X3FdmuK6PbznTRqwTEv3NXKVUzHrzzTd57733eOCBB2JyoFqbVq9ezcSJE7n99tttR1GWGG8Rq996jKKCfE797d06t3eQ/bDgf/n27cfI2fSN7SgnTBuJSqmItHr1am666SYGDRrElClTbMeJOT179uSuu+7iH//4BzNmzLAdR1mwad6r7NmwnC4XTqZ2k1a248Sc9qOuoVajFnzz2gPk7d9tO84J0UaiUioiZWRkMHz4cN577z0SEvTKmFB45JFHGDFiBDfffDNLly61HUeFWZ1mp5B51m9odcb5tqPEpITkNHpf/xjeAg8rX76HIk++7UjHTX/yKqUiypEjR0hOTqZ169Z89NFHtuPENLfbzbvvvkufPn0YP3482dnZekNLDKlonL6CvMMkptamWc/BNOs5OMyp4kudppn0vPJPrJhxNxv+9RJdLpxc4b9LJF6rqGcSlVIRIzc3l/PPP59LL71U5xgOk3r16vHvf/+bV199VRuIceDAtg189tAl/Lz6C9tR4kaTbgPoddWDtDsn+maK0kaiUioi5OTkMHToUBYuXMjYsWN1xocwatu2LcOHDwfgjTfe4JtvovdCe1WxPeuXs3TqZBJS0sho1cl2nLjSrNcQkmrVxVtYwKZ5r1KYHx3zO2sjUSll3bp16xgwYACrVq1i5syZXH755bYjxaWjR4/y0EMPcdZZZ/Hhhx/ajqOCxBjD1sUfsfylO0lr0IwBt79Aar3GtmPFpb0/rGbTv19j6dTJ5O3/xXacKmkjUSllVUFBAaNHj+bAgQNkZWUxbtw425HiVkpKCosWLaJjx46MGzeOSZMmkZcXHWc8VMX2fr+Kte89RcP2p3HGrc+RUreh7Uhxq1GH3vS5/jEO79rKf564OuK7/bWRqJSy4rvvvsPj8ZCYmMg777zD6tWrOfPMyLtwO96cdNJJLF68mN///vc8//zz9OrVC48n+qcXizder5eDP2UD0KBdT3rf8Dh9J/6VxNTalpOpJt0GMPDOGaQ1bM7X/7iPjXNfsR2pQtpIVEqF1bfffssll1xC165dee655wDo168fTZs2tZxMHZOcnMzf/vY3PvvsM2666SaSkopnkFi4cCFer9dyOlWZwsJC3n77bbp3786ip28ib99uRISm3c6MiWniYkXtxq0YcPsLdBh9A01LpkLcvn07mzZtspwskNXvGBG5Q0SMiOi5b6WOg4iMFJGNIpItIneV87qIyLMlr68RkV42ch5jjOHFF1/kjDPOoGfPnsybN48pU6Zw1VXRd7dfPBk8eDC33XYbAEuXLmXIkCG0b9+eBx54gA0bNsT0HejRVmM7d+7kj3/8I61bt/Zd03vqZXeTktHIZixVCZc7gXYjrqRuyw4APPbYY3Ts2JERI0bw2muvceDAAbsBsdhIFJGWwHBgm60MSkUjEXED04BRQGfgUhHpXGqzUUC7kuVG4IVw5SsqKmLjxo289957TJ069Vhm3nnnHQ4dOsSTTz7Jtm3bePzxx2nQoEG4Yqka6tOnD++88w6tW7fmkUceoVOnTrRp04b169cDcODAAY4ePWo5ZXBEeo0dOXKEpUuX8sILL7Bw4UIAPB4Pzz77LL1792b27NmsWbOG5qcN07OHUeTBBx/k/vvvJzs7m2uuuYaGDRtywQUX+F7/5Zdfwv6Hmc3BtJ8B7gR0tFyljk9fINsYsxlARN4FxgLfObYZC7xhin+iLBORDBE5yRjzc2U7NsaQm5tLQUGBb2nUqBEJCQnk5OSwY8cOfv31Vw4cOOD7/8SJE0lKSuKZZ55h+vTpbN261XezQ61atbjhhhtISUnh448/pk6dOjq0TZRKSEhgwoQJTJgwgZ07d/LBBx/w+eefk5mZCcATTzzBX//6VzIzM2nVqhUtW7akVatWPPzww4gIq1evZv/+/aSnp5OcnExSUhKpqam0aNECwHfdo9vtxmW/YROyGvN4PAH1JSI0bFjcmfb9999z4MCBgPpq2rQpo0ePBmD48OFs3LiRHTt2+BoLV199NYMHDyYzM5NffvmFunXrBvmjUOHSpEkTHnroIR588EG++uor5syZ47vUA6BXr14cPHiQNm3a+GpsxIgRvobkv//9b9LT00lLSyM5OZnk5GTq169PRkYGxhjy8/Nxu9243dWfo9tKI1FEzgd+Msas1l8YSh235sB2x+MdwOnV2KY5UOkvsG+++YZatWoFPLd582batGnDjBkzuPvuu8u8Z8KECTRu3JiMjAy6dOnCqFGj6NatGz179qRz586+H3Lp6enV/PJUpGvWrBmTJk1i0qRJvufOO+88EhISyM7OZtu2bSxcuJCDBw/yyCOPAPDkk0/yv//7vwH7ady4Mbt3F89pO378eGbPng3ARRddFKavpEIhqbH169eTnJwc8NyAAQNYtKh4Bo6xY8f6zsweM2LECF8jsUmTJjRr1oyTTz6ZHj160KNHD1q18s+5rA3E2CAi9OvXj379+vmeM8bwwAMPsHbtWrZs2cL27dtZsmQJABdccAGFhYWMGjWqzL7+8Ic/8NRTT3Ho0KET+v6QUJ26FJEsoLwr0e8F7gFGGGN+FZEtQG9jTE4F+7mR4lP5AF2B/4YgbjA1BMr9WiJEpOeD6MjYwRhTx8aBRWQ8cI4x5vqSx1cAfY0xtzi2mQM8boxZVPJ4AXCnMebrcvYXTTUWDd8bmjE4YqLGoqy+IDq+NzRjzVWrvkJ2JtEYM6y850WkG9AGOHYWsQXwjYj0NcbsKmc/04HpJe9daYzpHarMwRDpGSM9H0RPRouH3wG0dDxuAew8gW2A6KqxSM8HmjFYYqXGoqm+QDMGS6RnrG59hf3CD2PMWmNMY2NMpjEmk+JC61VeA1EpVa4VQDsRaSMiScAEYHapbWYDV5bcgdkP+LWqa6WUUj5aY0ph98YVpdQJMMYUisjvgE8AN/CKMWadiEwsef1FYC5wLpAN5ALX2MqrVLTRGlOqmPVGYsnZxOqaHqocQRTpGSM9H2jGKhlj5lL8S8r53IuOdQNMKv2+aoj0zz7S84FmDJZYrDH93INDM9ZctfKF7MYVpZRSSikVvawPRqWUUkoppSJP1DYSI3VKPxH5q4hsKJmm6QMRybCd6ZiqppmyTURaishCEVkvIutE5FbbmcojIm4RWSUi/7KdJVQitb5Aa+xERUt9gdaYbZFaY5FcXxCbNRaVjUSJ7Cn95gNdjTHdgU1A2dGHLZDqTTNlWyHwB2NMJ6AfMCkCMwLcCqyvcqsoFeH1BVpjJypa6gu0xmyLuBqLgvqCGKyxqGwk4p/SL+IuqDTGfGqMKSx5uIzisbMigW+aKWOMBzg2zVTEMMb8bIz5pmT9EMXfwM3tpgokIi2A84AZtrOEUMTWF2iNnahoqC/QGosEEVpjEV1fEJs1FnWNRHFM6Wc7SzVcC8yzHaJERVNIRSQRyQR6Al9ZjlLa3yn+4e61nCMkoqy+QGvshERwfYHWWKSJlBqLmvqC2Kkx60PglEeqMaVfeBMFqiyfMeajkm3upfjU89vhzFaJ8ibJjsi/YkWkNjATuM0Yc9B2nmNEZDTwizHmaxE523KcExbp9QVaY6EUqfUFWmPhFIU1FhX1BbFVYxHZSAzWlH7hzneMiFwFjAaGmsgZY6ja07TZJCKJFBfX28aYWbbzlDIAOF9EzgVSgHQRecsYc7nlXMcl0usLtMZCJcLrC7TGwiYKayzi6wtir8aiepxEEdkC9DbGRMwk2iIyEngaGGSM2WM7zzEikkDxBchDgZ8onnbqMmPMOqvBHKT4p+brwD5jzG2W41Sq5C+wO4wxoy1HCZlIrC/QGjtR0VRfoDVmUyTWWKTXF8RmjUXdNYlR4DmgDjBfRL4VkRerekM4lFyEfGyaqfXAPyOpuEoMAK4AhpR8dt+W/LWjlJPW2InR+lLVFXE1FgX1BTFYY1F9JlEppZRSSoWGnklUSimllFJlaCNRKaWUUkqVoY1EpZRSSilVhjYSlVJKKaVUGdpIVEoppZRSZWgjUSmllFJKlaGNRKWUUkopVYY2EuOEiPQRkTUikiIitURknYh0tZ1LqVig9aVU6IlIpohsEJHXS+rt/0QkzXauWKaDaccREXmU4rkaU4EdxpjHLUdSKmZofSkVWiKSCfwInGmMWSwirwDfGWOespssdmkjMY6ISBLF810eBc4wxhRZjqRUzND6Uiq0ShqJ/zHGtCp5PASYbIy5wGauWKbdzfGlPlCb4jk5UyxnUSrWaH0pFXqlz2zpma4Q0kZifJkO3A+8DTxhOYtSsUbrS6nQayUi/UvWLwUW2QwT6xJsB1DhISJXAoXGmP8VETewRESGGGM+s51NqWin9aVU2KwHrhKRl4DvgRcs54lpek2iUkoppSJeyTWJ/zLG6MgBYaLdzUoppZRSqgw9k6iUUkoppcrQM4lKKaWUUqoMbSQqpZRSSqkytJGolFJKKaXK0EaiUkoppZQqQxuJSimllFKqDG0kKqWUUkqpMv4/7oekkJoLrDgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 792x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3)\n",
    "axs[0].hist2d(x_list, v_list, bins=51, range=[[-4, 4], [-4, 4]],\n",
    "              cmap=\"Blues\", density=True)\n",
    "axs[0].set_xlim(-4, 4)\n",
    "axs[0].set_ylim(-4, 4)\n",
    "axs[0].set_xlabel(\"x\")\n",
    "axs[0].set_ylabel(\"p\")\n",
    "axs[0].set_xticks([-4, -2, 0, 2, 4])\n",
    "axs[0].set_yticks([-4, -2, 0, 2, 4])\n",
    "axs[0].set_aspect(1)\n",
    "\n",
    "position = np.linspace(-4, 4, 500)\n",
    "boltz_factor = np.exp(-0.5/kT*k_spring*position**2)\n",
    "Z = np.trapz(boltz_factor, position)\n",
    "boltz_factor /= Z\n",
    "axs[1].plot(position, boltz_factor, color=\"k\", linestyle=\"--\")\n",
    "axs[2].plot(position, boltz_factor, color=\"k\", linestyle=\"--\")\n",
    "\n",
    "axs[1].hist(x_list, density=True, bins=21, alpha=0.8)\n",
    "axs[1].set_xlim(-4, 4)\n",
    "axs[1].set_ylim(0, 0.5)\n",
    "axs[1].set_box_aspect(1)\n",
    "axs[1].set_xlabel(\"x\")\n",
    "axs[1].set_ylabel(\"probability\")\n",
    "axs[1].set_xticks([-4, -2, 0, 2, 4])\n",
    "\n",
    "axs[2].hist(v_list, density=True, bins=21, alpha=0.8)\n",
    "axs[2].set_xlim(-4, 4)\n",
    "axs[2].set_ylim(0, 0.5) \n",
    "axs[2].set_box_aspect(1)\n",
    "axs[2].set_xlabel(\"p\")\n",
    "axs[2].set_ylabel(\"probability\")\n",
    "axs[2].set_xticks([-4, -2, 0, 2, 4])\n",
    "\n",
    "fig.set_size_inches(11, 3)\n",
    "plt.savefig(\"md_ho_nh.pdf\", bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ovito",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
